﻿#include "RBF_Class.h"

using namespace std;

int main()
{
    //    interpolation_test2(); // use invert function to calculate polynomial coefficients

    ofstream outfile("F:\\VisualStudio_WorkSpace\\output\\output_Bessel\\Bessel_2lambda_Eigen.log");
    if (!outfile) {
        std::cout << "outfile open fail" << std::endl;
    }

    const double miu0 = 1.2566370614359172953850573533118e-6;
    const double epsilon0 = 8.85418781762039e-12;
    const double c_wave_speed = 1.0 / sqrt(miu0 * epsilon0);
    std::cout << "c_wave_speed = " << c_wave_speed << "\n";
    double freq = 30E+09;
    const double lambda = c_wave_speed / freq;
    //const double lambda = 0.333333333;
    int Nums_1D = 3, Max_Nums_1D = 10;
    int Start_i = 0, Max_i = 20, Max_j = 20, Max_k = 20;
    double c_new;
    vector<std::pair<vector<int>, vector<double>>> Nums1Dxyz_Lengthxyz_OptimalC_Error_list;
    double RBF_min_Average_error = 0.0;
    pair<double, double> RBF_LBF_Error;
    vector<pair<vector<int>, double>> Pre_i_RBF_Error_list, Pre_j_RBF_Error_list;

    for (int i = 0; i < Max_i; i++)
    {
        double length_x = (0.1 + i * 0.1) * lambda;
        for (int j = 0; j < Max_j; j++)
        {
            double length_y = (0.1 + j * 0.1) * lambda;
            for (int k = 12; k < Max_k; k++)
            {
                double length_z = (0.1 + k * 0.1) * lambda;
                outfile << "length_x = " << (0.1 + i * 0.1) << "   length_y = " << (0.1 + j * 0.1) << "   length_z = " << (0.1 + k * 0.1) << "\n";
                for (int Nums_1D_x = 4; Nums_1D_x <= Max_Nums_1D; Nums_1D_x++)
                {
                    for (int Nums_1D_y = 6; Nums_1D_y <= Max_Nums_1D; Nums_1D_y++)
                    {
                        for (int Nums_1D_z = 9; Nums_1D_z <= Max_Nums_1D; Nums_1D_z++)
                        {
	                        // ReSharper disable once StringLiteralTypo
	                        std::cout << "////// TartanGrids" << "   Nums_1D_x = " << Nums_1D_x << "   Nums_1D_y = " << Nums_1D_y << "   Nums_1D_z = " << Nums_1D_z << "//////" << "\n";

	                        try
	                        {
		                        RBF rbf;
		                        RBF_LBF_Error = rbf.New_test_scattered_RBF(length_x, length_y, length_z, lambda, TartanGrids,
		                                                                   Nums_1D_x, Nums_1D_y, Nums_1D_z, -1,
		                                                                   Nums1Dxyz_Lengthxyz_OptimalC_Error_list,
		                                                                   Bessel);
                            }
                            catch (const std::exception& e)
                            {
                                std::cerr << "Exception caught: " << e.what() << "\n";
                            }
                            catch (...)
                            {
                                std::cerr << "Unknown exception caught" << "\n";
                            }

                            if (RBF_LBF_Error.first < 0.02 || abs(RBF_LBF_Error.first - 0.02) < 1E-10)
                            {
                                outfile << "Nums_1D_x = " << Nums_1D_x << "   Nums_1D_y = " << Nums_1D_y << "   Nums_1D_z = " << Nums_1D_z << "\n";
                                double optimal_c = Nums1Dxyz_Lengthxyz_OptimalC_Error_list.back().second[3];
                                double RBF_error = Nums1Dxyz_Lengthxyz_OptimalC_Error_list.back().second[4];
                                double Lagrange_error = Nums1Dxyz_Lengthxyz_OptimalC_Error_list.back().second[5];
                                outfile << "optimal_c = " << optimal_c << "\n";
                                outfile << "RBF_error = " << RBF_error << "   Lagrange_error = " << Lagrange_error << "\n";
                                outfile << "\n====================================================================\n\n";
                                break;
                            }
                            if (Nums_1D_x == Max_Nums_1D && Nums_1D_y == Max_Nums_1D && Nums_1D_z == Max_Nums_1D && RBF_LBF_Error.first > 0.02)
                            {
                                outfile << "---------------------------------------------------\n";
                                outfile << "All interpolation points do not meet the 0.02 error\n";
                                outfile << "Nums_1D_x = " << Nums_1D_x << "   Nums_1D_y = " << Nums_1D_y << "   Nums_1D_z = " << Nums_1D_z << "\n";
                                double optimal_c = Nums1Dxyz_Lengthxyz_OptimalC_Error_list.back().second[3];
                                double RBF_error = Nums1Dxyz_Lengthxyz_OptimalC_Error_list.back().second[4];
                                double Lagrange_error = Nums1Dxyz_Lengthxyz_OptimalC_Error_list.back().second[5];
                                outfile << "optimal_c = " << optimal_c << "\n";
                                outfile << "RBF_error = " << RBF_error << "   Lagrange_error = " << Lagrange_error << "\n";
                                outfile << "---------------------------------------------------\n\n";
                            }
                        }
                        if (RBF_LBF_Error.first < 0.02 || abs(RBF_LBF_Error.first - 0.02) < 1E-10)
                        {
                            break;
                        }
                    }
                    if (RBF_LBF_Error.first < 0.02 || abs(RBF_LBF_Error.first - 0.02) < 1E-10)
                    {
                        break;
                    }
                }
            }
        }
    }
    std::cout << "Test finished" << "\n";
    return 0;
    //new_test_scattered_lagrange(); //use SVD function to calculate polynomial coefficients
}